﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Diplomova_prace
{
    class Integration_alg
    {
      //  Polynoms p = new Polynoms();
        Expression v = new Expression();
        Integrace_ret ret = new Integrace_ret();
        List<List<double>> weights = new List<List<double>>();
        List<List<double>> nodes = new List<List<double>>();
        public delegate void PostHandler();
        public event PostHandler post;
        bool incorrect = false;

        public Integration_alg()
        {
            v.incorrectInput += v_incorrectInput;
        }

        void v_incorrectInput()
        {
            if (!incorrect)
                post();
            incorrect = true;
        }
        public Integrace_ret MidpointRule(string function, double a, double b)
        {
            incorrect = false;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            string parser = v.RPN(function);
            double fs = v.ComputeRPN(parser, (a + b) / 2);
            double h = b - a;
            fx.Add(fs);
            x.Add((a + b) / 2);
            ret.fxProperty = fx;
            ret.xProperty = x;
            ret.widthProperty = h;
            ret.returnProperty = h * fs;
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            ret.stepProperty = h;
            return ret;

        }
        public Integrace_ret ComMidpointRule(string function, double a, double b, int numOfPart)
        {
            incorrect = false;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            string parser = v.RPN(function);
            double h = (b - a) / numOfPart;
            ret.numOfPartProperty = numOfPart;
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            double functionalValue = 0;
            double valueOfPoint;
            for (int i = 0; i < numOfPart; i++)
            {
                valueOfPoint = v.ComputeRPN(parser, a+h/2.0);
                x.Add(a + h / 2.0);
                fx.Add(valueOfPoint);
                functionalValue += valueOfPoint;
                
                a = a + h;
         
            }
            ret.widthProperty = h;
            ret.returnProperty = h * functionalValue;
            ret.xProperty = x;
            ret.fxProperty = fx;
            ret.stepProperty = h;
            return ret;
        }
        public Integrace_ret TrapezodialRule(string function, double a, double b)
        {
            incorrect = false;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            string parser = v.RPN(function);
            double fa = v.ComputeRPN(parser, a) ;
            double fb = v.ComputeRPN(parser, b) ;
            double fs = (fa + fb);
            double h = (b - a);
            fx.Add(fa);
            fx.Add(fb);
            x.Add(a);
            x.Add(b);
            ret.fxProperty = fx;
            ret.xProperty = x;
            ret.widthProperty = h / 2.0;
            ret.returnProperty = (h / 2.0) * fs;
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            ret.stepProperty = h;
            return ret;
        }
        public Integrace_ret ComTrapezodialRule(string function, double a, double b, int numOfPart)
        {
            incorrect = false;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            string parser = v.RPN(function);
            double h = (b - a) / numOfPart;
            ret.numOfPartProperty = numOfPart;
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            double functionalValue = 0;
            double valueOfPoint;
            for (int i=0; i <= numOfPart; i++)
            {
                valueOfPoint = v.ComputeRPN(parser, a);
                x.Add(a);

                if (i == 0 || i == numOfPart)
                {
                    functionalValue += valueOfPoint;
                    fx.Add(valueOfPoint);
                }
                else
                {
                    functionalValue += 2 * valueOfPoint;
                    fx.Add(2*valueOfPoint);
                }
                a = a + h;
               
            }
            ret.widthProperty = h /2.0;
            ret.returnProperty = h/2.0 * functionalValue;
            ret.xProperty = x;
            ret.fxProperty = fx;
            ret.stepProperty = h;
            return ret;
            //return vys;
            //vysledek = (h / 2) * funkcniHodnoty;
            //return vysledek;
        }
        public Integrace_ret SimpsonRule(string function, double a, double b)
        {
            //dva dily standartne
            incorrect = false;
            string parser = v.RPN(function);
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            double h = (b - a);
            double fa = 0; //funkcni hodnoty
            double fs = 0;
            double fb = 0;
            fa = v.ComputeRPN(parser, a);
            fs = v.ComputeRPN(parser, (a + b) / 2);
            fb = v.ComputeRPN(parser, b);
            fx.Add(fa);
            fx.Add(fs);
            fx.Add(fb);
            x.Add(a);
            x.Add((a + b)/ 2);
            x.Add(b);
            ret.fxProperty = fx;
            ret.xProperty = x;
            ret.widthProperty = h / 6.0;
            ret.returnProperty = (h / 6.0) * (fa + 4 * fs + fb);
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            ret.stepProperty = h;
            return ret;
         // h/3, h=(b-a)/2, casteji se uvadi h/6 coz odpovida h/3, h=(b-a)/2 -> h/6
        }
        public Integrace_ret ComSimpsonRule(string function, double a, double b, int numOfPart)
        {
            
            //pocet dilu musi byt sudy -> lichy pocet uzlu - kvuli neprekryvani intervalu 
            incorrect = false;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            string parser = v.RPN(function);
            if (numOfPart % 2 != 0)
                throw new Exception("Počet dílů musí být sudý");
            ret.numOfPartProperty = numOfPart;
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            ret.widthProperty = (b - a) / (3 * numOfPart);
            double h = (b - a) / (numOfPart); 
      //      double funkcniHodnoty = 0; //funkcni hodnoty
            double result = 0;
            b = a + 2*h;
            double fa = 0; //funkcni hodnoty
            double fs = 0;
            double fb = 0;
            
            for (int i = 0; i < numOfPart/2; i++)
            {
                // h/3, h=(b-a)/2, casteji se uvadi h/6 coz odpovida h/3, h=(b-a)/2 -> h/6
                fa = v.ComputeRPN(parser, a);
                fs = v.ComputeRPN(parser, (a + h));
                fb = v.ComputeRPN(parser, b);
                fx.Add(fa);
                fx.Add(fs);
                x.Add(a);
                x.Add((a + h));
                
                result += ((b - a) / 6.0) * (fa + 4 * fs + fb);
                a = b;
                b = b+2*h;
            }
            x.Add(b-2*h);
            fx.Add(v.ComputeRPN(parser, b - 2 * h));
            ret.returnProperty = result;
            ret.xProperty = x;
            ret.fxProperty = fx;
            ret.stepProperty = h;
            return ret;
        }
        public Integrace_ret BoolRule(string function, double a, double b)
        {
            incorrect = false;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            string parser = v.RPN(function);
            double h = (b - a);
            double fa =  v.ComputeRPN(parser, a);
            double fQuarter = v.ComputeRPN(parser, (3 * a + b) / 4);
            double fHalf =  v.ComputeRPN(parser, (a + b) / 2);
            double fThreeQuarters =  v.ComputeRPN(parser, (a + 3 * b) / 4);
            double fb =  v.ComputeRPN(parser, b);
            fx.Add(fa);
            fx.Add(fQuarter);
            fx.Add(fHalf);
            fx.Add(fThreeQuarters);
            fx.Add(fb);
            x.Add(a);
            x.Add((3 * a + b) / 4);
            x.Add((a + b) / 2);
            x.Add((a + 3 * b) / 4);
            x.Add(b);
            ret.widthProperty = h / 90;
            ret.returnProperty = (h / 90) * (7 * fa + 32 * fQuarter + 12 * fHalf + 32 * fThreeQuarters + 7 * fb);
            ret.xProperty = x;
            ret.fxProperty = fx;
            ret.upperBoundaryProperty = b;
            ret.lowerBoundaryProperty = a;
            ret.stepProperty = h;
            return ret;
        }
        public Integrace_ret SlozeneBooleovoPravidlo(string function, double a, double b, int numOfPart)
        {//skripta chyba ve vzorci - > sirka 2*h/45 !!!
                       //pocet dilu musi byt sudy -> lichy pocet uzlu - kvuli neprekryvani intervalu 
            incorrect = false;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            string parser = v.RPN(function);
            if (numOfPart % 4 != 0)
                throw new Exception("Počet dílů musí být dělitelný čtyřmi");
            double result = 0;
            double h = (b - a) / (numOfPart);
          
            ret.numOfPartProperty = numOfPart;
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            ret.widthProperty = 2*(b - a) / (45 * numOfPart);
            b = a + 4 * h;
            double fa;
            double fQuarter;
            double fHalf;
            double fThreeQuarters;
            double fb;
            for (int i = 0; i < numOfPart / 4 ; i++)
            {
                fa = v.ComputeRPN(parser, a);
                fQuarter = v.ComputeRPN(parser, (3 * a + b) / 4);
                fHalf = v.ComputeRPN(parser, (a + b) / 2);
                fThreeQuarters = v.ComputeRPN(parser, (a + 3 * b) / 4);
                fb = v.ComputeRPN(parser, b);
                fx.Add(fa);
                fx.Add(fQuarter);
                fx.Add(fHalf);
                fx.Add(fThreeQuarters);
         //       fx.Add(fb);
                x.Add(a);
                x.Add((3 * a + b) / 4);
                x.Add((a + b) / 2);
                x.Add((a + 3 * b) / 4);
            //    x.Add(b);
                result += ((b - a) / 90) * (7 * fa + 32 * fQuarter + 12 * fHalf + 32 * fThreeQuarters + 7 * fb);
                a = b;
                b = a + 4 * h;
               
            }
            x.Add(b-4*h);
            fx.Add(v.ComputeRPN(parser, b - 4 * h));
            ret.returnProperty = result;
            ret.xProperty = x;
            ret.fxProperty = fx;
            ret.stepProperty = h;
            return ret;
        }
        public Integrace_ret ThreeEightsRule(string function, double a, double b)
        {
            incorrect = false;
            string parser = v.RPN(function);
            double h = (b - a);
       //     double vysledek = 0;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            double fa = v.ComputeRPN(parser, a);
            double fThird =  v.ComputeRPN(parser, (2 * a + b) / 3.0);
            double fTwoThirds =  v.ComputeRPN(parser, (a + 2 * b) / 3.0);
            double fb = v.ComputeRPN(parser, b);
            fx.Add(fa);
            fx.Add(fThird);
            fx.Add(fTwoThirds);
            fx.Add(fb);
            x.Add(a);
            x.Add((2 * a + b) / 3.0);
            x.Add((a + 2 * b) / 3.0);
            x.Add(b);
            ret.fxProperty = fx;
            ret.xProperty = x;
            ret.widthProperty = h / 8.0;
            ret.returnProperty = (h / 8.0) * (fa + 3 * fThird + 3 * fTwoThirds + fb);
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            ret.stepProperty = h;
            return ret;
          //  vysledek = (h / 8) * fx; // h/3, h=(b-a)/2, casteji se uvadi h/6 coz odpovida h/3, h=(b-a)/2 -> h/6
        }
        public Integrace_ret ComThreeEights(string function, double a, double b, int numOfPart)
        {
            incorrect = false;
            List<double> fx = new List<double>();
            List<double> x = new List<double>();
            string parser = v.RPN(function);
            double h = (b - a) / (numOfPart);
            ret.numOfPartProperty = numOfPart;
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            ret.widthProperty = (3*(b - a)) / (8 * numOfPart);
            double result = 0;
            b = a + 3 * h;
            double fa;
            double faa;
            double fbb;
            double fb;
            for (int i = 0; i < numOfPart / 3; i++)
            {
               fa = v.ComputeRPN(parser, a);
               faa =  v.ComputeRPN(parser, (2 * a + b) / 3 );
               fbb =  v.ComputeRPN(parser, (a + 2 * b)  / 3);
               fb = v.ComputeRPN(parser, b);
               fx.Add(fa);
               fx.Add(faa);
               fx.Add(fbb);
             
               x.Add(a);
               x.Add((2 * a + b) / 3.0);
               x.Add((a + 2 * b) / 3.0);
            
               result += ((b-a)/8.0) * (fa + 3 * faa + 3 * fbb + fb);
               a = b;
               b += 3 * h;
            }
            x.Add(b - 3 * h);
            fx.Add(v.ComputeRPN(parser, b - 3 * h));
            ret.returnProperty = result;
            ret.xProperty = x;
            ret.fxProperty = fx;
            ret.stepProperty = h;
            return ret;
        }
        public Integrace_ret RombergIntegration(string function, double a, double b, double accuracy)
        {
            incorrect = false;
            double[,] trapezodial = new double[20, 20];
            for (int i = 0; i < trapezodial.GetLength(0); i++)
            {
                trapezodial[i, 0] = ComTrapezodialRule(function, a, b, (int)Math.Pow(2, i)).returnProperty;
                for (int j = 1; j <= i; j++)
                {
                    trapezodial[i, j] = ((Math.Pow(4, j) * trapezodial[i, j - 1] - trapezodial[i - 1, j - 1]) / (Math.Pow(4, j) - 1));
                }
                if ( i > 0)
                    if (Math.Abs(trapezodial[i, i] - trapezodial[i, i - 1]) < accuracy)
                    {
                        ret.RombergProperty = Copy(trapezodial, i);
                        return ret;
                    }
            }
            ret.RombergProperty = Copy(trapezodial, trapezodial.GetLength(0) - 1);
            return ret;
        }
        public Integrace_ret GaussQuadrature(string function, double a, double b, int numOfPart, int countOfNode)
        {
            // skripta chyba ve vzorci, vnitrni suma jde 0 do n-1 !!!
            string ba = v.RPN(function);
            double ab = v.ComputeRPN(v.RPN(function), -4);
                incorrect = false;
            string parser = v.RPN(function);
            GaussWeightsAndNodes();
            double h = (b - a) / numOfPart;
            double fx = 0;
            double actualfx;
            List<double> currfx = new List<double>();
            double tmp;
            for (int i = 0; i < numOfPart; i++)
            {
                for (int j = 0; j < countOfNode; j++)
                {
                    tmp = v.ComputeRPN(parser, 0.5 * nodes[countOfNode - 2][j] * h + a + h * (i + 0.5));
                    actualfx = weights[countOfNode - 2][j] * v.ComputeRPN(parser, 0.5 * nodes[countOfNode - 2][j] * h + a + h * (i + 0.5)); 
                    fx += actualfx;
                    currfx.Add(tmp);

                }
            }
            ret.widthProperty = h/2.0;
            ret.lowerBoundaryProperty = a;
            ret.upperBoundaryProperty = b;
            ret.GaussPointProperty = nodes[countOfNode - 2];
            ret.GaussWeightProperty = nodes[countOfNode - 2];
            ret.fxProperty = currfx;
            ret.numOfPartProperty = numOfPart;
            ret.returnProperty = fx * h / 2.0;
            return ret;

        } 
        public  void GaussWeightsAndNodes ()
        {
            weights.RemoveRange(0, weights.Count);
            nodes.RemoveRange(0, nodes.Count);
        List<double> weightsGauss2 = new List<double>() {1,1 };
        List<double> weightsGauss3 = new List<double>() { 8.0 / 9.0, 5.0 / 9.0, 5.0 / 9.0 };
        List<double> weightsGauss4 = new List<double>() { 1.0 / 36.0 * (18 + Math.Sqrt(30)), 1.0 / 36.0 * (18 + Math.Sqrt(30)), 1.0 / 36.0 * (18 - Math.Sqrt(30)), 1.0 / 36.0 * (18 - Math.Sqrt(30)) };
        List<double> weightsGauss5 = new List<double>() { 128.0 / 225.0, 1.0 / 900.0 * (322 + 13 * Math.Sqrt(70)), 1.0 / 900.0 * (322 + 13 * Math.Sqrt(70)), 1.0 / 900.0 * (322 - 13 * Math.Sqrt(70)), 1.0 / 900.0 * (322 - 13 * Math.Sqrt(70)) };
        List<double> xi2 = new List<double>() { (1 / (Math.Sqrt(3))), (1 / (-Math.Sqrt(3))) };
        List<double> xi3 = new List<double>() { 0, 1.0 / 5.0 * (Math.Sqrt(15)), -1.0 / 5.0 * (Math.Sqrt(15)) };
        List<double> xi4 = new List<double>() { 1.0 / 35.0 * Math.Sqrt(525 - 70 * Math.Sqrt(30)), -1.0 / 35.0 * Math.Sqrt(525 - 70 * Math.Sqrt(30)), -1.0 / 35.0 * Math.Sqrt(525 + 70 * Math.Sqrt(30)), 1.0 / 35.0 * Math.Sqrt(525 + 70 * Math.Sqrt(30)) };
        List<double> xi5 = new List<double>() { 0, 1.0 / 21.0 * Math.Sqrt(245 - 14 * Math.Sqrt(70)), -1.0 / 21.0 * Math.Sqrt(245 - 14 * Math.Sqrt(70)), 1.0 / 21.0 * Math.Sqrt(245 + 14 * Math.Sqrt(70)), -1.0 / 21.0 * Math.Sqrt(245 + 14 * Math.Sqrt(70)) };
        weights.Add(weightsGauss2);
        weights.Add(weightsGauss3);
        weights.Add(weightsGauss4);
        weights.Add(weightsGauss5);
        nodes.Add(xi2);
        nodes.Add(xi3);
        nodes.Add(xi4);
        nodes.Add(xi5);
        }
        public double[,] Copy(double[,] maticeA, int rozmer)
        {
            double[,] maticeB = new double[rozmer+1,rozmer+1];
            for (int i = 0; i <= rozmer; i++)
            {
                for (int j = 0; j <= rozmer; j++)
                {
                    maticeB[i, j] = maticeA[i, j];
                }
            }
            return maticeB;
        }

    }
    class Integrace_ret
    {
        double width;
        double result;
        double lowerBoundary;
        double horniMez;
        double numOfPart;
        double step;
        List<double> x;
        List<double> fx;
        List<double> GaussVahy;
        List<double> GaussBody;
        double[,] Romberg = new double[20, 20];

        public double stepProperty
        {
            get { return step; }
            set { step = value; }
        }
        public double lowerBoundaryProperty
        {
            get { return lowerBoundary; }
            set { lowerBoundary = value; }
        }
        public double upperBoundaryProperty
        {
            get { return horniMez; }
            set { horniMez = value; }
        }
        public double widthProperty
        {
            get{ return width;}
            set { width = value; }
        }
        public double returnProperty
        {
            get { return result; }
            set { result = value; }
        }
        public List<double> xProperty
        {
            get { return x; }
            set { x = value; }
        }
        public List<double> fxProperty
        {
            get { return fx; }
            set { fx = value; }
        }
        public double numOfPartProperty
        {
            get { return numOfPart; }
            set { numOfPart = value; }
        }
        public List<double> GaussWeightProperty
        {
            get { return GaussVahy; }
            set { GaussVahy = value; }
        }
        public List<double> GaussPointProperty
        {
            get { return GaussBody; }
            set { GaussBody = value; }
        }
        public double[,] RombergProperty
        {
            get { return Romberg; }
            set { Romberg = value; }
        }
    }
}
